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New  York  University 

New  York,  New  York  10012 


ABSTRACT 

The  compressible  Rayle-gh-Taylcr  instability  of  a  supersonic  accelerated  con- 
tact discontinuity  between  two  gases  is  studied  by  numerically  solving  the 
two-dimensional  Euler  equations.  The  computed  solutions  exhibit  a  compli- 
cated set  of  nonlinear  waves  comprised  of  spike  and  bubble  bow  shocks,  ter- 
minal shocks  within  the  spike  and  bubble,  Kelvin-Helmholtz  rollup  of  the 
spike  tip,  and  contact  surface  waves.  The  spike  appears  to  attain  a  finite 
growth  of  aspect  ratio  approximately  equal  to  2.  The  propagation  of  a  super- 
sonic slab  jet  is  also  studied  numerically,  in  order  to  compare  and  contrast  the 
jet  wave  structure  with  that  of  the  supersonic  accelerated  surface. 


1.  This  work  was  supported  by  tlie  Applied  Mathematical  Sciences  subprogram  of  the  Office  of  Energy 
Research,  U.S.  Department  of  Energy  under  Contraa  DE-AC02-76ER03077. 


1.  Introduction 

The  classical  instabilities  of  hydrodynamics  exhibit  interesting  new  features  in  the  context  of 
supersonic  gas  dynamics.  In  the  present  investigation,  we  study  the  compressible  Rayleigh-Taylor 
instability  of  a  supersonic  accelerated  contact  discontinuity  between  two  gases.   We  also  study  the 
propagation  of  a  supersonic  slab  jet  in  order  to  compare  and  contrast  the  jet  wave  structure  with  that 
of  the  supersonic  accelerated  surface.   These  supersonic  instabilities  and  flows  lead  to  the  formation 
of  complicated  sets  of  shock  waves  and  to  the  appearance  of  Kelvin-Helmholtz  rollup  and  contact 
surface  waves. 

In  this  paper,  we  examine  the  evolution  into  the  deeply  nonlinear  regime  of  a  single  mode  in 
the  Rayleigh-Taylor  instability  (a  single  finger).   In  a  subsequent  paper,  we  plan  to  examine  the 
statistics  of  several  modes  (many  fingers).   The  focus  of  the  present  investigation  is  to  validate  the 
use  of  front-tracking  on  coarse  grids  for  supersonic  surface  instabilities,  as  a  prelude  to  investigating 
the  statistical  or  chaotic  surface  instability  regime. 

These  processes  arise  in  a  range  of  physical  applications:  the  accelerated  siurface  problem 
models  important  features  of  the  Rayleigh-Taylor  and  Meshkov-Richtmeyer  instabilities,  with 
applications  to  inertial  laser  fusion,  instabilities  of  layers  in  stars,  and  the  instability  of  laser 
accelerated  foils;  the  supersonic  jet  results  can  be  applied  to  laboratory,  aircraft,  and  astrophysical 

jets. 

The  computations  are  performed  using  a  front-tracking  code  [1]  which  has  been  validated  on  a 
wide  range  of  problems  for  which  independent  comparison  solutions  or  experimental  data  exist 
including  propagation  of  an  expanding  or  contracting  circularly  symmetric  shock  wave,  steady-state 
supersonic  flow  past  a  wedge,  Kelvin-Helmholtz  vortex  rollup  [1],  regular  reflection  of  a  shock 
wave,  and  Mach  reflection  of  a  shock  wave  [2].   An  incompressible  code  using  the  same  front- 
tracking  method  has  been  validated  on  two-phase  flow  in  oil  reservoirs  [3-5]  and  on  the 
incompressible  Rayleigh-Taylor  instability  [6]  (for  which  high  quality  numerical  comparison  solutions 
[7,8]  are  known). 
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The  code  employs  a  finite  difference  method  together  with  the  tracking  of  selected  waves 
(contact  discontinuities  in  this  investigation)  to  solve  the  two  dimensional  Euler  equations  of 
compressible  gas  dynamics  in  conservation  form.  The  Euler  equations  for  a  compressible,  inviscid, 
polytropic  gas  can  be  cast  in  the  form  of  a  general  hyperbolic  system  of  nonlinear  conservation  laws 
(1.1)        w^  +  Vf(w)  =  0 
by  setting 


(1.2) 


■(;) 


and  f(w)  = 


m 
mm 


p  is  the  mass  density,  m  is  the  momentimi  density,  E  is  the  total  energy  density,  and  p  is  the 
thermodynamic  pressure  specified  by  the  poly  tropic  equation  of  state 


(1.3) 


(^-1) 


|mr 
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These  equations  express  the  conservation  of  mjiss,  momentum,  and  energy. 

The  front-tracking  method  substantially  improves  the  space-time  grid  resolution  of  our 
computations.   The  one  dimensional  grids  tracking  the  fronts  evolve  dynamically  and  are 
automatically  refined,  leading  to  enhanced  resolution  of  the  solution  in  the  entire  computational 
region.   Consider  a  series  of  runs  with  increasingly  finer  meshes  covering  the  computational  region. 
As  the  mesh  is  refined  with  the  number  of  tracked  waves  held  constant,  more  of  the  difficulty  of 
computing  the  solution  is  transferred  to  the  interior  scheme,  and  the  advantages  of  front  tracking 
may  be  diminished.   However,  if  the  interface  is  unstable  ,  there  is  an  increase  in  detail  in  the 
convolutions  of  the  interface  as  the  mesh  is  refined  and  a  corresponding  increase  in  the 
computational  difficulty  of  resolving  the  interface,  so  that  the  front-tracking  method  retains  its 
advantages.   This  enhanced  resolution  is  essential  in  resolving  many  fingers  within  a  reasonable 
amount  of  computational  time,  since  to  resolve  many  fingers  accurately  on  a  practical  computational 
grid,  a  computational  method  must  be  able  to  resolve  a  single  finger  accurately  on  a  coar'.e  grid. 

We  demonstrate  the  increased  resolution  of  the  front-tracking  code  by  presenting  supersonic 
jet  results  on  40x60  and  80x120  grids  and  comparing  them  with  the  more  detailed  results  [9]  of 


Norman,  Smarr,  and  Winkler  (NSW)  on  a  160x300  grid.  NSW  use  the  method  of  LeBlanc  interface 
tracking,  which  is  a  fractional  volume  interface-tracking  method,  to  obtain  accurate  resolution  of  the 
contact  boundary.  Using  the  "surface"  front-tracking  method  [1],  we  are  able  to  get  good  results 
(including  contact  shape  and  wave  structure)  with  a  jet  beam  only  2.5  grid  blocks  across.  This 
capability  appears  to  be  unique  to  the  front-tracking  method. 

The  accelerated  surface  and  supersonic  jet  runs  have  been  validated  under  mesh  refinement  by 
factors  up  to  four  in  each  linear  dimension  for  the  entire  range  of  Mach  numbers  and  density  ratios 
presented  here.   In  addition,  the  wave  structure  of  the  computed  solutions  and  the  morphology  of  the 
contact  discontinuity  were  analyzed  and  compared  with  the  wave  structure  of  the  results  of  NSW  and 
with  expected  scaling  behavior  of  the  Rayleigh-Taylor  instability  as  a  function  of  density  ratio  for  the 
single  finger. 

2.  Supersonic  accelerated  surfaces 

A  compressible  gas  interface  which  is  accelerated  by  a  shock  (the  Meshkov-Richtmeyer 
instability  [10-12])  or  by  a  strong  gravitational  field  is  Rayleigh-Taylor  unstable  when  the  gases  are 
accelerated  toward  each  other.   The  important  features  of  this  instability  can  be  modelled  by 
imparting  an  initial  kinetic  energy  to  the  contact  discontinuity,  which  subsequently  is  allowed  to 
advect  freely.   We  imagine  that  the  source  of  the  acceleration  is  a  shock  wave  or  a  gravitational  force 
acting  during  a  short  time,  after  which  the  acceleration  is  turned  off.   To  analyze  this  instability,  we 
imagine  that  a  contact  surface  between  a  light  (above)  and  heavy  gas  (below)  is  perturbed  as  a  sine 
wave  with  one  mode,  with  a  velocity  field  given  initially  by 
it(TAcos(fcT)exp(CTf)cosh(a(y->'^^„)) 


(2.1)        u  = 


V  = 


asinh(a(Tio->'ft^„)) 
crA  sin(fcj:)exp(cTr)sinh(o(y -yi,rf„,)) 


sinh(a(Tio->'wo)) 

where  «  is  the  x  velocity,  v  is  the  y  velocity,  ti  =  Aexp(CTr)sin(fcc)  is  the  y  position  of  the  surface  in 
the  linearized  theory,  tiq  =  Ti(r=0)  ,  a  is  the  growth  rate  of  the  surface,  k  is  the  wave  number  of  the 

perturbation,  a  =  ^k?-^<j^/c^  is  the  wave  number  in  the  y  direction  ,  and  ^-^^    is  the  top  (bottom)  of 


the  computational  domain.  We  assume  that  the  problem  is  periodic  in  x  with  reflecting  boundaries  at 

y  =  vwo-  • 

The  problem  can  be  parametrized  in  dimensionless  units  by  the  initial  Mach  number  of  the  tip 
of  the  spike  with  respect  to  the  heavy  gas  and  by  the  initial  density  ratio  D  -  Pt/Pg  {b  denotes  the 
gas  below  the  contact,  a  the  gas  above).   The  dimensional  scales  are  set  by  the  initial  ambient 
pressure,  perturbation  wavelength,  and  initial  ambient  density  of  the  heavy  gas.   V/e  choose  the 
ambient  density  of  the  heavy  gas  so  that  the  sound  speed  in  the  ambient  heavy  gas  equals  1.  The 
sound  speed  in  the  ambient  heavy  gas  and  the  perturbation  wavelength  define  a  dimensionless  time 
scale.   The  polytropic  gas  constant  -y  was  set  equal  to  1.4.   Our  computations  sample  the  parameter 
space  defined  by  density  ratios  equal  to  2,  8,  and  100,  and  by  Mach  numbers  equal  to  2.8  (£>  =  2),  3 
(D  =  8),  and  2.5  {D  =  100).  These  Mach  numbers  indicate  the  upper  limits  for  the  given  density 
Trtios  beyond  which  our  solution  method  fails,  presumably  due  to  oscillations  in  the  interior  scheme 
produced  by  strong  waves  and/or  to  an  initialization  which  is  too  far  removed  from  the  Meshkov 
shock-induced  initialization.   For  validation  purposes,  we  also  give  results  for  Mach  number  0.5  and 
density  ratio  2. 

The  problem  can  be  initialized  with  a  linearized  solution  of  the  compressible  Euler  equations 
(wave  equation)  or  the  velocity  field  perturbation  described  by  Eq.  (2.1)  above,  depending  on  the 
physical  situation  to  be  modelled.   For  Mach  numbers  less  than  or  equal  to  0.5,  the  two  initializations 
give  very  similar  results,  and  agree  with  linear  theory  for  small  amplitudes.   The  subsequent 
nonlinear  d.-^velopment  of  the  surface  is  largely  independent  of  the  initialization  for  a  given  Mach 
number  of  the  spike.   Figure  1  compares  the  interfaces  for  Mach  number  0.5  and  density  ratio  2  at 
f  =  5  for  the  two  initializations.   For  Mach  numbers  greater  than  or  equal  to  1,  the  Euler  equations 
do  not  linearize  since  0  ^  |a|,|v|  s  A/  ,  and  the  velocity-field  perturbation  described  by  Eq.  (2.1) 
should  be  used.   For  most  of  our  nms  we  used  a  grid  of  40x80  points.   Our  results  were  validated 
under  grid  refinement  on  a  mesh  of  80x160  points.   For  M  =  2.8  and  D  =  2  comparisons  could  be 
made  among  20x40,  40x80,  and  80x160  grids.   Tlie  initialization,  Euler  equations,  and  boundary 
conditions  respect  Ici't-right  symmetry  about  the  vertical  center-line.   Tiiis  symmetry,  however,  is  not 


imposed  in  the  computations.  Thus  the  observed  left-right  symmstry  of  the  computed  solutions 
serves  as  a  check  on  the  computations. 

An  interesting  set  of  wave  structures  emerges  from  this  study.  We  have  employed  plots  of  the 
contact  discontinmty  along  with  pressure  and  density  contour  plots  to  exhibit  these  wave  structures. 
In  addition,  cross  section  plots  of  density  and  pressure  parallel  to  the  y  axis  along  the  center  and  edge 
of  the  computational  region  were  made  to  highlight  the  relative  strengths  of  the  shock  and  rarefaction 

waves. 

Typically  as  the  spike  advances  into  the  light  gas,  it  develops  the  characteristic  Rayleigh-Taylor 
vortex  rollup.   The  advancing  spike  generates  a  bow  wave  or  a  bow  shock  in  the  surrounding  light 
gas,  depending  on  whether  the  flow  is  subsonic  or  supersonic  in  the  light  gas.   As  the  bubble  at 
either  side  of  the  computational  domain  advances  into  the  heavy  gas,  bow  waves  or  shocks  are 
generated  in  the  heavy  gas.  The  bow  shocks  in  front  of  the  bubbles  and  the  spike  end  in  zero 
strength.  The  bow  shocks,  after  transporting  kinetic  energy  away  from  the  contact,  effectively 
decouple  from  the  subsequent  development  of  the  contact  instability.   (In  the  case  of  reflecting 
boundeiry  conditions  at  the  top  and  bottom  of  the  computational  region,  the  bow  shocks  produce 
reflected  waves  which  interact  with  the  contact  at  late  times  and  influence  its  shape  slightly.) 
Eventually  the  bubble  bow  shocks  intersect  and  cross,  and  the  spike  bow  shock  interacts  with  its 
periodic  neighbors. 

The  bow  shocks  are  actually  N  waves  (a  rarefaction  with  a  shock  on  each  side)  with  the 
trailing  shock  of  zero  strength  (see  Figure  2),  which  result  from  the  impulsive  initialization.  The 
bubble  and  spike  behave  analogously  to  a  piston  which  is  instantaneously  given  a  finite  velocity  and 
then  slowed  down  by  hitting  ambient  gas,  producing  a  shock  followed  by  a  rarefaction  of  almost 
equal  strength.   N  waves  have  a  decay  rate  proportional  to  1/ vr.  They  are  physical  transients  or 
next  to  leading  order  in  the  large  time  asymptotic  expansion.  The  leading  order  in  the  asymptotic 
expansion  is  given  by  the  net  difference  of  total  shock  strength  minus  total  rarefaction  strength  in  the 
N  wave.   This  difference  is  left  after  all  wave  cancellation,  and  does  not  decay  at  all.   The  net 
difference  in  the  computed  N  waves  is  very  small  compared  to  the  strength  of  the  bow  shock  or  the 


rarefaction  wave. 

Just  inside  of  the  advancing  spike  a  "terminal"  shock  wave  is  formed.   The  contact  is  advancing 
more  slowly  than  the  heavy  gas  inside  of  the  spike.   A  compression  or  shock  wave,  preceded  by  a 
rarefaction  wave,  is  formed  as  the  advancing  heavy  gas  is  slowed  down  to  the  contact  velocity.   The 
same  phenomenon  can  be  seen  in  supersonic  jets  [9],  and  was  described  theoretically  by  Blanford  and 
Rees  [13]  in  their  original  paper  on  astrophysical  jets.  Once  the  terminal  wave  has  formed  inside  the 
spike,  it  propagates  in  the  negative  y  direction  until  it  interacts  with  the  reflected  conjoined  bubble 
bow  wave  (for  reflecting  boimdary  conditions  at  the  lower  y  boundary).   If  through-flow  boundary 
conditions  are  imposed  at  the  lower  y  boundary,  first  the  conjoined  bubble  bow  wave  and  then  the 
terminal  wave  propagate  out  of  the  computational  domain.   A  similcU'  terminal  wave  is  formed  in  the 
light  gas  inside  of  the  advancing  bubble. 

This  shock  preceded  by  a  rarefaction  wave  pattern  can  be  dearly  seen  in  the  center-line  and 
edge  plots  of  density  in  both  the  accelerated  surface  mns  (Figures  2,  3,  and  6)  and  the  jet  runs 
(Figures  9  and  10).   Note  that  while  the  jet  terminal  shock  propagates  with  the  head  of  the  jet  beam, 
the  accelerated  surface  terminal  shocks  are  physical  transients  which  decouple  from  the  late  evolution 
of  the  contact  instability.   However,  the  terminal  shocks  do  play  an  important  role  in  the  early 
development  of  the  instability. 

Figures  2-4  portray  the  detailed  evolution  of  the  Mach  2.8  density  ratio  2  accelerated  surface. 
The  flow  is  initially  supersonic  in  both  gases.   At  time  t  =  0.4  ,  the  bow  shocks  in  the  lower  gas 
have  interacted  to  form  a  single  shock,  while  the  spike  bow  shock  has  interacted  and  joined  with  its 
periodic  neighbors.   At  late  times  {t  =  0.8),  the  spike  exhibits  highly  developed  rollup,  and  the 
contact  shape  indicates  the  presence  of  small-scale  surface  instabilities.   The  center-line  and  edge-line 
density  plots  indicate  the  location  of  the  contact  and  the  relative  strengths  of  the  bow  shocks  and 
terminal  shocks.   Note  that  the  pattern  of  waves  (moving  along  the  positive  y  direction)  for  the 
center-line  density  cross  section  is  repeated  (mutatis  mutandis,  since  the  waves  are  wcciker  in  the  light 
gas)  in  reverse  for  the  edge  cross  section. 

Figure  5  compares  the  results  for  the  Mach  2.8  density  ratio  2  accelerated  surface  obtained  on 
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20x40,  40x80,  and  80x160  grids. 

Similar  structure  can  be  seen  in  the  density  ratio  8  Mach  2.5  and  density  ratio  100  Mach  2.5 
runs.  Results  are  presented  in  Figures  6  and  7.  Note  that  for  a  given  aspect  ratio  of  spike  height  to 
perturbation  wavelength,  the  finger  becomes  thinner  and  the  spike  develops  less  rollup  with 
increasing  density  ratio,  as  is  true  for  the  incompressible  Rayleigh-Taylor  instability. 

The  compressible  Rayleigh-Taylor  results  differ  from  the  incompressible  case  chiefly  in  the 
formation  of  the  terminal  compression  waves  and  in  the  fact  that  the  spike  exhibits  less  rollup.   The 
accelerated  surface  problem  differs  from  the  gravitational  instability  in  that  the  spike  appears  to 
attain  a  finite  growth  of  aspect  ratio  approximately  equal  to  2  for  our  range  of  parameters.  The 
accelerated  surface  is  imparted  an  initial  kinetic  energy  which  is  converted  first  into  the  initial  shock 
waves  and  then  into  "turbulent"  eddies  at  the  rontact  boundary  and  radiated  sound  waves  on  a  time 
scale  of  about  2.   Thus  the  spike  velocity  rapidly  decreases  and  the  bubble  appears  to  come  almost  to 
rest.   Here  we  do  not  intend  to  make  any  statement  about  the  motion  of  the  spike  or  bubble  when 
the  flow  becomes  essentially  incompressible  at  very  low  Mach  numbers  «  0.1  due  to  the  fact  that 
there  are  no  dissipative  forces  for  steady  flow  in  an  ideal  incompressible  fluid  (D'Alembert's 
paradox).   Plots  of  spike  and  bubble  velocity  versus  time  for  Mach  number  2.8,  density  ratio  2  are 
presented  in  Figure  8.   In  the  incompressible  (gravitational)  Rayleigh-Taylor  case  [7,8],  the  spike 
achieves  a  terminal  free-fall  acceleration,  while  the  bubble  achieves  a  terminal  velocity,  due  to  the 
constant  addition  of  gravitational  energy  to  the  system.   We  believe  that  the  finite  growth  of  the 
spike  in  the  accelerated  surface  problem  is  reproduced  in  the  Meshkov-Richtmeyer  instability,  since 
the  Meshkov  unstable  contact  also  is  imparted  a  finite  kinetic  energy  by  the  incident  shock,  and  then 
advects  freely. 

3.  Supersonic  jets 

Supersonic  jets  were  £imong  the  first  compressible  phenomena  to  be  investigated.   A 
resurgence  of  interest  has  been  sparked  by  the  observation  of  astrophysical  jets  emanating  from  the 
cores  of  active  galaxies  and  by  the  subsequent  success  of  theoretical  and  computational  analyses.   We 
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have  made  a  series  of  slab  jet  nins  and  have  compared  our  40x60  and  80zl20  grid  results  to  the 
160x300  grid  results  of  NSW  [9].   The  jet  beam  in  our  computations  is  2.5  i;, .id  blocks  wide  for  the 
40x60  grid  and  5  grid  blocks  wide  for  the  80x120  grid,  while  the  beam  is  20  grid  blocks  wide  in  the 
computation  of  NSW.   (NSW  use  70  ratioed  zones  on  either  side  of  the  jet  beam.) 

Almost  all  of  the  supersonic  jets  of  NSW  are  axisymmetric.  Although  our  compressible  code 
at  present  lacks  cylindrical  symmetry,  we  plan  to  install  cylindrical  symmetry  and  will  report  on  the 
results  elsewhere. 

The  jet  was  initialized  by  injecting  gas  at  a  specified  Mach  number  into  an  ambient  gas  at  equal 
pressure.   The  boundary  conditions  are  through-flow  (coupling  to  an  ambient  reservoir  specified  by 
tlie  initialization),  which  allow  waves  to  propagate  out  of  the  computational  domain.   The  problem  is 
parametrized  by  the  Mach  number  of  th;?  jr.t  witVi  '•espect  to  the  jet  gas  and  the  density  ratio  of  jet  to 
ambient  gas.  We  choose  the  initial  density  of  the  jet  £as  so  that  the  sound  speed  equals  1  in  the 
undisturbed  jet  gas.   -y  was  set  equal  to  5/3,   The  problem  is  independent  of  length  scale,  so  the 
results  apply  to  jets  from  laboratory  to  astrophysical  scales. 

We  have  studied  in  detail  the  evolution  of  a  Mach  3,  density  ratio  10  slab  jet  (see  Figure  9). 
At  present  we  are  unable  to  reach  Mach  12  (the  published  slab  jet  results)  with  our  interior  solver. 
Nonetheless,  the  similarities  of  our  Mach  3  slab  jet  with  the  Mach  12  slab  jet  of  Ref.  [9]  are  striking. 

Figure  9  indicates  the  presence  of  a  bow  wave  (the  flow  is  subsonic  in  the  ambient  gas)  and  of 
a  terminal  shock  near  the  head  of  the  jet  beam,  preceded  by  a  rarefaction  wave.   The  more  detailed 
computations  of  NSW  indicate  that  the  terminal  shock  is  in  fact  a  Mach  disc  with  associated  incident 
and  reflected  shocks  and  slip  lines.   This  terminal  shock  system  may  explain  the  observed  hot  spots 
terminating  astrophysical  jets  [9].   The  contact  shape  displays  the  large  scale  Kelvin- Helmholtz  rollup 
of  this  jet,  and  the  development  of  two-dimensional  pinch  waves. 

Figures  9  and  10  compare  the  development  of  the  jet  at  t  =  0.4  for  an  80x120  and  a  40x60 
grid.   The  density  contours  indicate  the  computed  wave  structure,  and  the  center  line  plots  indicate 
the  relative  strength  of  the  rarefaction  wave  and  the  terminal  sliock.   Tne  jet  is  2.5  grid  blocks  v/ide 
in  Figure  10  and  5  grid  blocks  wide  in  Figure  9. 
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The  fact  that  we  get  reasonable  results  with  a  beam  2.5  grid  blocks  across  ill"-strates  one  of  the 
advantages  of  the  front-tracking  method.   By  placing  additional  degrees  of  freedom  around  the 
tracked  contact,  the  method  is  able  to  resolve  the  solution  globally  with  fewer  degrees  of  freedom 
than  required  by  conventional  finite  difference  methods,  since  in  two  dimensions  if  the  computational 
grid  has  0{N^)  points,  the  tracked  front  has  0{N)  points.  The  importance  of  this  feature  of  the 
method  will  become  apparent  when  multiple  fingers  are  considered. 

4.  A  descriptioQ  of  the  front-tracking  method 

Front-tracking  is  an  adaptive  computational  method  for  solving  a  hyperbolic  system  of  non- 
linear conservation  laws.  In  two-dimensional  problems,  a  moving  one-dimensional  griu,  called  the 
front,  is  fitted  to  and  tracks  selected  waves  in  the  solution.  These  waves  can  be  ^^.arp  discontinuities 
which  exist  as  mathematical  solutions  of  idealized  physical  equations  (e.g.  the  Euler  equations)  or 
waves  for  which  physical  quantities  change  rapidly  but  smoothly  over  a  fraction  of  a  mesh  length 
(e.g.  chemical  reaction  fronts).   For  compressible  fluid  dynamics,  these  waves  include  shock  waves, 
contact  discontinuities,  materifd  interfaces,  phase  boundaries,  slip  lines,  and  chemical  reaction  fronts. 

The  front  divides  the  computational  grid  into  topologically  connected  interior  regions  called 
components.  The  solution  is  computed  by  first  propagating  the  front  and  then  the  solution  in  each 
component. 

The  front  is  advanced  in  two  steps.   First  the  Rankine-Hugoniot  equations  are  used  to 
propagate  the  front  normally  by  solving  a  non-local  Riemann  problem.  Then  tangential  waves  are 
propagated  along  the  front  using  a  one- dimensional  Lax-Wendroff  method.   At  points  (called  nodes) 
where  discontinuity  curves  intersect,  the  propagation  is  defined  by  the  solution  of  shock  polar 
equations,  as  a  first  approximation  to  solving  a  two-dimensional  Riemann  problem.   The  propagation 
of  the  front  is  described  more  thoroughly  in  Ref.  [1],  while  front  tracking  and  two-dimensional 
Riemann  problems  are  discussed  in  Ref.  [2]. 

The  interior  regions  between  fronts  are  treated  as  initijil/boundary  problems  and  the  solutions 
in  these  regions  are  computed  using  an  operator  split  Lax-Wendroff  finite-difference  method.   The 
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front  and  interior  schemes  are  coupled  in  a  strip  of  width  one  mesh  spacing  on  either  side  of  the 
front. 
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Figure  Captions 

Fig.  1.   Comparison  of  M  =  0.5,  D  =  2  computed  contact  fronts  at  t  =  5  which  evolved  from  the 

velocity  field  and  compressible  wave  equation  initializations.   20x40  grid. 

Fig.  2.  M  =  2.8,  D  =  2  density  contours  and  density  cross  sections  at  t  =  0.2.  80x160  grid. 

Fig.  3.   M  =  2.8,  D  =  2  density  contours  and  density  cross  sections  at  t  =  0.4.   80x160  grid. 

Fig.  4.   M  =  2.8,  D  =  2  contact  at  t  =  0.8.   80x160  grid. 

Fig.  5.   Comparison  of  M  =  2.8,  D  =  2  contacts  at  t  =  0.8.   (a)  20x40  and  80x160  grids,   (b)  40x80 

and  80x160  grids. 

Fig.  6.   M  =  3,  D  =  8  density  contours  and  density  cross  sections  at  t  =  0.4.   80x160  grid. 

Fig.  7.   M  =  2.5,  D  =  100  contact  at  t  =  0.3.   80x160  grid. 

Fig.  8.   Spike  and  bubble  velocities  vs.  t  for  M  =  2.8,  D  =  2.   40x160  grid  with  through-flow 

boundary  conditions.  -^    '-• 

Fig.  9.   Density  contours  and  density  cross  section  for  M  =  3,  D  =  10  jet  at  t  =  0.4.   80x120  grid. 

Fig.  10.   Density  contours  and  density  cross  section  for  M  =  3,  D  =  10  jet  at  t  -  0.4.   40x60  grid. 
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